Skrypty geoprzetwarzania

Analizy geoprzestrzenne

Author

Krzysztof Dyba

# install.packages("remotes")
remotes::install_github("kadyb/rgugik")
library("terra")
library("rgugik")

Główny Urząd Geodezji i Kartografii jest istotnym źródłem danych przestrzennych dla Polski. Dane można przeglądać i pobrać z Geoportalu lub wykorzystując różne usługi. W otwartych zbiorach danych znajdziemy m. in.:

Wyszukiwanie i pobieranie wymienionych zbiorów danych umożliwia pakiet rgugik.

Ortofotomapa

Ortofotomapa to rastrowe, ortogonalne i kartometryczne przedstawienie powierzchni terenu powstałe w wyniku cyfrowego przetwarzania zdjęć lotniczych lub satelitarnych. Podczas ortorektyfikacji usuwane są zniekształcenia geometryczne wynikające z rzeźby terenu przy użyciu cyfrowych modeli wysokości. Ortofotomapa posiada georeferencje, co pozwala na określenie współrzędnych geograficznych dla każdej komórki obrazu.

Cechy ortofotomapy:

  • Rozdzielczość przestrzenna – związana z rozmiarem najmniejszego obiektu, który może zostać wykryty przez czujnik i jest określana przez rozmiar komórki obrazu (piksel). Im mniejsza komórka, tym więcej szczegółów reprezentuje. Zbyt duży rozmiar oznacza, że poszczególne obiekty na zdjęciu przestają być rozpoznawalne.
  • Kompozycja – obrazy analogowe są przedstawione w odcieniach szarości, natomiast obrazy cyfrowe mogą składać się z naturalnych kolorów (RGB) lub bliskiej podczerwieni (NIR).

Wyszukiwanie

Do wyszukania arkuszy ortofotomapy służy funkcja ortho_request(), która jako argument przyjmuje geometrię. Jako przykładową geometrię możemy wykorzystać punkt. W tym celu należy stworzyć macierz, w której wierszy reprezentują punkty (w naszym przypadku tylko jeden), a kolumny współrzędne X i Y. Następnie należy dokonać konwersji do obiektu wektorowego używając funkcji vect() i definiując odpowiedni układ współrzędnych.

punkt = cbind(16.92, 52.41)
punkt = vect(punkt, crs = "EPSG:4326")
dane = ortho_request(punkt)

Możemy wyświetlić część otrzymanej ramki danych lub alternatywnie przeglądać całość używając funkcji View().

# wyświetl 10 pierwszych wierszy i 6 pierwszych kolumn
dane[1:10, 1:6]
            sheetID year resolution composition  sensor        CRS
1  N-33-130-D-d-1-2 2023       0.25         CIR Digital    PL-1992
2   6.177.11.05.3.3 2023       0.05         CIR Digital PL-2000:S6
3        6.177.11.2 2023       0.50         RGB Digital PL-2000:S6
4   6.177.11.05.3.3 2023       0.05         RGB Digital PL-2000:S6
5     6.177.11.05.3 2023       0.10         CIR Digital PL-2000:S6
6     6.177.11.05.3 2023       0.10         RGB Digital PL-2000:S6
7       6.177.11.05 2023       0.20         CIR Digital PL-2000:S6
8       6.177.11.05 2023       0.20         RGB Digital PL-2000:S6
9        6.177.11.2 2023       0.50         CIR Digital PL-2000:S6
10 N-33-130-D-d-1-2 2021       0.25         CIR Digital    PL-1992

Standardowo dane możemy filtrować z uwzględnieniem zadanych parametrów.

dane[dane$year > 2016, ]
dane[dane$composition == "CIR", ]

I sortować, np. według aktualności.

# kolejność malejąca (najnowsze dane)
dane[order(-dane$year), ]

Pobieranie

Jako przykład pobierzmy dwie kompozycje tego samego obszaru wykonane w naturalnych barwach i z kanałem bliskiej podczerwieni z 2021 r. (ID: 75107_1047046_N-33-130-D-d-1-2 i 75106_1052150_N-33-130-D-d-1-2).

id = c("75107_1047046_N-33-130-D-d-1-2", "75106_1052150_N-33-130-D-d-1-2")
dane_sel = dane[dane$filename %in% id, ]

Po selekcji potrzebnych danych, można je pobrać wykorzystując funkcję tile_download(). Możliwe jest również wskazanie katalogu, do którego powinny zostać pobrane obrazy (argument outdir).

Zazwyczaj warto zwiększyć domyślną wartość przekroczenia czasu połączenia (timeout) z domyślnych 60 sekund w przypadku dużych plików lub wolnego połączenia.

options(timeout = 600)
tile_download(dane_sel, outdir = "dane")

Do wylistowania pobranych plików służy funkcja list.files(). Należy wskazać jakie pliki chcemy wczytać (pattern = "\\.tif$") i zapobiegawczo zwrócić pełne ścieżki do plików (full.names = TRUE).

pliki = list.files("dane", pattern = "\\.tif$", full.names = TRUE)
pliki
[1] "dane/75106_1052150_N-33-130-D-d-1-2.tif"
[2] "dane/75107_1047046_N-33-130-D-d-1-2.tif"

W ostatnim kroku możemy kolejno wczytać rastry i je wyświetlić.

# kompozycja w naturalnych barwach
r1 = rast(pliki[1])
plot(r1)

# kompozycja z bliską podczerwienią
r2 = rast(pliki[2])
plot(r2)

Cyfrowe modele wysokościowe

Cyfrowe modele wysokościowe to modele opisujące powierzchnię terenu. Powstają w wyniku obróbki zdjęć lotniczych, skanowania laserowego (LiDAR), pomiarów geodezyjnych czy interferometrii radarowej (InSAR). CMW są jednym z kluczowych zbiorów danych w systemach informacji geograficznej i stanowią podstawę wielu środowiskowych analiz przestrzennych. Ponadto są źródłem produktów pochodnych, takich jak cieniowanie, nachylenie czy szorstkość terenu.

CMW to ogólna nazwa grupy modeli o różnych cechach, uwzględniając:

  • Numeryczny model terenu (Digital Terrain Model) – reprezentacja pozbawiona jakichkolwiek obiektów nad powierzchnią terenu, takich jak budynki czy drzewa.
  • Numeryczny model pokrycia terenu (Digital Surface Model) – reprezentacja terenu wraz z jego pokryciem.

https://commons.wikimedia.org/w/index.php?title=File:DTM_DSM.svg

Cechy CMW:

  • Format – można wyróżnić trzy główne struktury: GRID (regularna siatka punktów / komórek), TIN (nieregularna topologiczna sieć trójkątów) oraz linie konturowe (dane wektorowe). Obecnie najczęściej używanym formatem jest GRID.
  • Dokładność – związana jest z pionowym błędem pomiaru.
  • Rozdzielczość przestrzenna – związana jest z rozmiarem najmniejszego obiektu, który może zostać wykryty przez czujnik i jest określana przez rozmiar komórki obrazu (piksel). Im większa komórka, tym bardziej uogólniona forma terenu.

Meteoryt Morasko

Naszym obszarem analiz jest rezerwat przyrody Meteoryt Morasko położony w północnej części Poznania. Został on utworzony w 1976 roku w celu ochrony obszaru kraterów uderzeniowych, które według badaczy powstały w wyniku upadku meteorytu Morasko około 5 tysięcy lat temu. Ponadto ochroną objęty jest las grądowy z rzadkimi gatunkami roślin i ptaków.

https://naukawpolsce.pl/aktualnosci/news%2C402631%2Csto-lat-temu-odkryto-pierwszy-fragment-meteorytu-morasko.html

(Fot. PAP © 2012 / Jakub Kaczmarczyk)

Więcej informacji znajdziesz tutaj:

Centroid (środek geometryczny) rezerwatu Morasko znajduje się na 16,895° długości geograficznej (X) i 52,489° szerokości geograficznej (Y).

wspolrzedne = matrix(c(16.895, 52.489), ncol = 2)
centroid = vect(wspolrzedne, type = "points", crs = "EPSG:4326")
centroid
 class       : SpatVector 
 geometry    : points 
 dimensions  : 1, 0  (geometries, attributes)
 extent      : 16.895, 16.895, 52.489, 52.489  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 

Dokonajmy również konwersji ze współrzędnych geograficznych na układ metryczny EPSG:2180.

centroid = project(centroid, "EPSG:2180")

W kolejnym kroku stwórzmy przybliżoną strefę, która obejmie obszar rezerwatu.

bufor = buffer(centroid, width = 400)

Stworzyliśmy bufor o szerokości 400 m, a teraz przygotujmy prostą wizualizację.

plot(bufor, main = "Bufor rezerwatu Morasko")
plot(centroid, col = "blue", add = TRUE)

Następnie możemy wyszukać dostępne dane wysokościowe dla tego obszaru za pomocą funkcji DEM_request() (jest ona analogiczna do funkcji ortho_request()). Jako argument musimy wskazać nasz utworzony poligon.

dane = DEM_request(centroid)

Oczywiście pozyskane dane możemy sprawdzić wywołując obiekt dane lub przejrzeć je interaktywnie używając funkcji View().

# wyświetl 10 pierwszych wierszy i 6 pierwszych kolumn
dane[1:10, 1:6]
                sheetID year    product              format         source
1           6.179.11.14 2021        DSM ARC/INFO ASCII GRID Laser scanning
2           6.179.11.14 2021        DTM ARC/INFO ASCII GRID Laser scanning
3      N-33-130-D-b-1-1 2021        DTM ARC/INFO ASCII GRID   Aerial photo
4      N-33-130-D-b-1-1 2020        DTM ARC/INFO ASCII GRID   Aerial photo
5  N-33-130-D-b-1-1-4-1 2019 PointCloud                 LAS Laser scanning
6      N-33-130-D-b-1-1 2012        DSM      ASCII XYZ GRID Laser scanning
7      N-33-130-D-b-1-1 2012        DSM ARC/INFO ASCII GRID Laser scanning
8      N-33-130-D-b-1-1 2012        DTM      ASCII XYZ GRID Laser scanning
9  N-33-130-D-b-1-1-4-1 2012 PointCloud                 LAS Laser scanning
10     N-33-130-D-b-1-1 2012        DTM ARC/INFO ASCII GRID Laser scanning
   resolution
1      0.50 m
2      0.50 m
3      5.00 m
4      5.00 m
5     12 p/m2
6      0.50 m
7      0.50 m
8      1.00 m
9     12 p/m2
10     1.00 m

Jak możemy zauważyć powyższe metadane opisują produkty o różnych formatach, aktualności, rozdzielczości oraz dokładności. Do naszej analizy potrzebujemy numerycznego modelu terenu (DTM) i numerycznego modelu pokrycia terenu (DSM) w formacie “ARC/INFO ASCII GRID”. Dokonajmy selekcji danych, tworząc dwie ramki danych i następnie łącząc je ze sobą przy pomocy funkcji rbind().

DTM_sel = dane[dane$format == "ARC/INFO ASCII GRID" &
               dane$product == "DTM" &
               dane$year == 2019, ]
DSM_sel = dane[dane$format == "ARC/INFO ASCII GRID" &
               dane$product == "DSM" &
               dane$year == 2019, ]

# połączenie powyższych ramek danych
dane_sel = rbind(DTM_sel, DSM_sel)
dane_sel[, 1:6]
            sheetID year product              format         source resolution
19 N-33-130-D-b-1-1 2019     DTM ARC/INFO ASCII GRID Laser scanning     1.00 m
15 N-33-130-D-b-1-1 2019     DSM ARC/INFO ASCII GRID Laser scanning     0.50 m

Wykorzystajmy funkcję tile_download() do pobrania tych dwóch produktów.

options(timeout = 600)
tile_download(dane_sel, outdir = "dane")

Przetwarzanie

Po pobraniu danych możemy je wczytać do sesji.

DTM = rast("dane/73044_917579_N-33-130-D-b-1-1.asc")
DSM = rast("dane/73043_917495_N-33-130-D-b-1-1.asc")

Możemy również nadać warstwom odpowiednie nazwy oraz, co ważniejsze, przypisać im układy współrzędnych zdefiniowane w obiekcie dane_sel (atrybut CRS).

# nadanie nazw warstwom
names(DTM) = "DTM"
names(DSM) = "DSM"

# ustawienie układu współrzędnych
crs(DTM) = crs(DSM) = "EPSG:2180"

Podczas pobierania prawdopodobnie zauważyłeś, że rastry różnią się czterokrotnie rozmiarem. Wynika to z różnicy ich rozdzielczości przestrzennej. W takiej sytuacji należy ujednolicić je do jednakowej rozdzielczości, aby móc je połączyć (nałożyć na siebie). Znacznie lepiej jest użyć niższej rozdzielczości niż ją sztucznie zwiększać, ponieważ nie możemy uzyskać więcej informacji, a przetwarzanie będzie szybsze. W tym celu użyjmy funkcji resample() i zapiszmy wynik na dysku określając ścieżkę w argumencie filename.

DSM = resample(DSM, DTM, method = "near", filename = "dane/DSM_1.tif")

Teraz oba modele mają te same wymiary (liczbę wierszy i kolumn) oraz rozdzielczość przestrzenną. Możemy więc połączyć je w jeden obiekt o nazwie DEM.

DEM = c(DTM, DSM)
nlyr(DEM)
[1] 2

Używając funkcji crop() możemy ograniczyć zasięg przestrzenny (bounding box) do zasięgu przestrzennego rezerwatu Morasko. W celu wykluczenia z analizy wartości poza obszarem poligonu, należy zastosować funkcję mask(). Można tę czynność również wykonać w jednym kroku, używają funkcji crop() z argumentem mask = TRUE.

DEM_crop = crop(DEM, bufor)
DEM_mask = mask(DEM_crop, bufor)

Zauważ, że domyślnie kolory zielony reprezentuje wysokie wartości, a kolor pomarańczowy niskie. Do wizualizacji danych topograficznych warto odwrócić paletę kolorów, tj. terrain.colors(99, alpha = NULL).

par(mfrow = c(1, 2)) # wyświetl 2 rastry obok siebie
plot(DEM_crop$DTM, main = "Docięcie", col = terrain.colors(99, alpha = NULL))
plot(bufor, add = TRUE)
plot(DEM_mask$DTM, main = "Maskowanie", col = terrain.colors(99, alpha = NULL))
plot(bufor, add = TRUE)

W pierwszej ćwiartce okręgu widzimy pięć mniejszych okręgów. Są to kratery powstałe po uderzeniu meteorytu Morasko. Największy znaleziony fragment waży 272 kg i jest największym meteorytem znalezionym w Polsce. Kolekcję można zobaczyć w Muzeum Ziemi WNGIG w Poznaniu.

Obliczmy szerokość krateru na podstawie poprzecznego profilu terenu. Przyjmijmy, że punkt A ma współrzędne (357280, 515980), a punkt B (357122, 515760).

punkty = matrix(c(357280, 515980,
                  357122, 515760),
                ncol = 2, byrow = TRUE)
linia = vect(punkty, type = "lines", crs = "EPSG:2180")
linia
 class       : SpatVector 
 geometry    : lines 
 dimensions  : 1, 0  (geometries, attributes)
 extent      : 357122, 357280, 515760, 515980  (xmin, xmax, ymin, ymax)
 coord. ref. : ETRF2000-PL / CS92 (EPSG:2180) 
plot(DEM_mask$DTM, main = "DTM [m]", col = terrain.colors(99, alpha = NULL))
plot(linia, col = "red", lwd = 2, add = TRUE)
text("A", x = 357320, y = 515980, cex = 0.8)
text("B", x = 357100, y = 515760, cex = 0.8)

W kolejnym kroku pozyskamy wartości wysokości dla wyznaczonego profilu za pomocą funkcji extract().

profil = extract(DEM, linia)
profil = profil$DTM
plot(profil, type = "l", xlab = "Indeks komórki", ylab = "Wysokość n.p.m. [m]")

Na profilu widoczny jest szum wynikający z metody pozyskania danych. Można to wygładzić w prosty sposób wykorzystując funkcję loess().

profil = loess(profil ~ seq_along(profil), span = 0.1)
profil = profil$fitted
plot(profil, type = "l", xlab = "Indeks komórki", ylab = "Wysokość n.p.m. [m]")

Powyższe ryciny na osi X zawierają indeks komórki. Lepszym rozwiązaniem będzie przedstawienie tej osi jako odległości od punktu startowego (A). Aby to wykonać musimy obliczyć średnią odległość między komórkami, czyli podzielić długość linii przez liczbę komórek.

odleglosc = perim(linia) / length(profil)
odleglosc
[1] 0.7146646

Następnie stwórzmy wektor odległości kolejnych punktów od punktu startowego wykorzystując sumę skumulowaną cumsum().

odleglosc = cumsum(rep(odleglosc, length(profil)))
odleglosc[1:5] # odległość pierwszych 5 punktów
[1] 0.7146646 1.4293293 2.1439939 2.8586585 3.5733232

Po tej operacji, możemy zaprezentować finalną wersję wykresu z zaznaczonym zasięgiem największego krateru uderzeniowego wynoszącą około 90 m.

plot(profil ~ odleglosc, type = "l", xlab = "Odległość [m]",
     ylab = "Wysokość n.p.m. [m]", main = "Numeryczny model terenu")
abline(v = c(50, 140), col = "blue")

W ostatnim kroku sprawdźmy wysokość obiektów (na tym obszarze są to drzewa). W tym celu należy odjąć NMT od NMPT. Różnica nazywana jest znormalizowanym NMPT, ponieważ przyjmuje wysokość terenu jako odniesienie.

nDSM = DEM$DSM - DEM$DTM
nDSM
class       : SpatRaster 
dimensions  : 2379, 2188, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 355732.5, 357920.5, 514649.5, 517028.5  (xmin, xmax, ymin, ymax)
coord. ref. : ETRF2000-PL / CS92 (EPSG:2180) 
source(s)   : memory
varname     : 73044_917579_N-33-130-D-b-1-1 
name        :       DSM 
min value   : -4.159996 
max value   : 67.769989 
# nadpisz wartości poniżej 0
nDSM[nDSM < 0] = 0

Powyższą operację odejmowania można wykonać również używając funkcji lapp().

plot(nDSM, main = "Wysokość drzew [m]",
     col = hcl.colors(9, palette = "Greens", rev = TRUE))

Zadanie

10. Pobierz minimum dwa sąsiadujące ze sobą kafelki ortofotomapy z tej samej serii i połącz je:

  1. do jednego pliku .tiff używając funkcji merge(),
  2. do jednego wirtualnego pliku .vrt używając funkcji vrt().

Sprawdź zajmowaną ilość miejsca przez te pliki na dysku wykorzystując funkcję file.size() (wynik zwracany jest w bajtach). Sprawdź również zawartość pliku .vrt (czym on jest w rzeczywistości?). Następnie, zmniejsz rozdzielczość mozaiki do 10 m i zapisz wynik. Jak zmieniła się jakość w stosunku do oryginalnego zdjęcia?

11. Wykonaj analizę ukształtowania terenu wybranego obszaru, w której uwzględnisz:

  • wizualizacje NMT i NMPT,
  • statystyki opisowe wysokości terenu (wartość minimalna, maksymalna, średnia oraz odchylenie standardowe),
  • profil wysokościowy,
  • wizualizacje wysokości obiektów (znormalizowany NMPT).

Pamiętaj, że:

  • Rastry powinny posiadać taką samą rozdzielczość przestrzenną oraz układ współrzędnych.
  • Różnica między NMPT a NMT nie może być ujemna.